<!DOCTYPE html>
<html lang="zh-cn">
<head>
    <title>数论算法</title>
    <meta charset="utf-8" />
    <link rel="stylesheet" type="text/css" href="../../css/note.css" />
</head>
<body>

<p class="algorithm">
  <b>n 次方根取整</b>
  求 `|__ root n a __|`, 其中整数 `a ge 0`, `n ge 2`.
  使用牛顿迭代法 `x = ((n-1)x + a//x^(n-1))//n`, 其中除法均向下取整.
</p>

<div class="playground" value="{ n: 2, a: 55555 }">
<p>说明: `235^2 = 55225`, `236^2 = 55696`, 因此 `|__ sqrt 55555 __| = 235`</p>
<script type="text">
module.exports = function floorRoot (str) {
  const { n, a } = Playground.parse(str)
  if (a < 0) return NaN
  let x = 1
  do x = Math.floor(
    ((n-1)*x + Math.floor(a/x**(n-1))) / n
  )
  while (x**n > a)
  return x
}
</script>
</div>

<p class="algorithm">
  <b>快速幂 (反复平方法)</b>
  求 `a^n`, 先将指数 `n` 写为二进制: `n = (n_k cdots n_1 n_0)_2`
  `= sum_(i=0)^k n_i 2^i`,
  于是
  <span class="formula">
    `a^n = prod_("where "n_i = 1) a^(2^i)`.
  </span>
</p>

<div class="playground" value="{ a: 5, n: 14, mod: Infinity }">
<p>
  在算法每次迭代的开头, 变量 `a = 5, 5^2, 5^4, 5^8`,
  对应到 `14` 的二进制表示 `0*1 + 1*2 + 1*4 + 1*8` 中系数为 `1` 的项,
  结果为 `5^14 = 5^2 * 5^4 * 5^8`.
</p>
<script type="text" id="pow">
function pow (a, n, mod = Infinity) {
  if (n < 0) return NaN
  let ret = 1
  while (n) {
    if (n & 1) ret = (ret * a) % mod
    a = (a * a) % mod
    n >>= 1
  }
  return ret
}

function demo (str) {
  let { a, n, mod } = Playground.parse(str)
  return pow(a, n, mod)
}

demo.pow = pow
module.exports = demo
</script>
</div>

<p class="algorithm">
  <b>辗转相除法</b> 设 a, b 为正整数, 求最大公约数 d,
  同时求出系数 x, y 满足 a x + b y = d.
</p>

<div class="playground" value="{ a: 75, b: 32 }">
<p>
75 = (1, 0)<br>
32 = (0, 1)<br>
11 = (1, -2)<br>
10 = (-2, 5)<br>
1 = (3, -7)
</p>
<script type="text" id="gcd">
function gcd (a, b) {
  if (b === 0) return a
  return gcd(b, a % b)
}

function gcdExtended (a, b) {
  if (b === 0) return { d: a, x: 1, y: 0 }
  const { d, x, y } = gcdExtended(b, a % b)
  // d = x * b + y * (a % b)
  //   = x * b + y * (a - b * (a // b))
  //   = y * a + (x - y * (a // b)) * b
  return { d, x: y, y: x - y * Math.floor(a / b) }
}

function demo (str) {
  const { a, b } = Playground.parse(str)
  return JSON.stringify(gcdExtended(a, b))
}

demo.gcd = gcd
demo.gcdExtended = gcdExtended
module.exports = demo
</script>
</div>

<p class="algorithm">
  <b>求模 `n` 乘法逆元</b>
  若 `a x -= 1 (mod n)`, 则称 `x` 为 `a` 模 `n` 的乘法逆元.<br>
  若 `(a, n) = 1`, 则存在 `x, y` 使得 `a x + n y = 1`, 即
  `a x -= 1 (mod n)`. 因此应用 gcdExtended 就可获得乘法逆.
</p>

<div class="playground" value="{ a: 42, n: 55 }">
<script type="text" id="inv">
const { gcdExtended } = require('gcd')
function inv (a, n) {
  let { x } = gcdExtended(a, n)
  return x < 0 ? x % n + n : x % n // 转换为最小非负剩余
}

function demo (str) {
  let { a, n } = Playground.parse(str)
  return inv(a, n)
}
demo.inv = inv
module.exports = demo
</script>
</div>

<p class="algorithm">
  <b>线性求逆</b>
  求 `1` 到 `p-1` 的各整数模素数 `p` 的乘法逆元.<br>
  设 `0 lt a lt p`, `x = p mod a`, `y = |__p // a__|`.
  则 `0 lt x lt p`, 说明 `x` 模 `p` 有逆元.
  计算知 `x + y a = p -= 0 (mod p)`,
  `x -= -y a (mod p)`,
  `x a^-1 -= -y (mod p)`,
  <span class="formula">
    `a^-1 = x^-1 (p-y) (mod p)`.
  </span>
</p>

<div class="playground" value="{ p: 17 }">
<script type="text">
module.exports = function linearInv (str) {
  let { p } = Playground.parse(str)
  const inv = [NaN, 1] // 0 没有逆; 1 的逆是 1
  for (let a = 2; a < p; ++a) {
    const [x, y] = [p % a, Math.floor(p / a)]
    inv[a] = inv[x] * (p-y) % p
  }
  return inv.slice(1)
}
</script>
</div>

<p class="algorithm">
  <b>前缀积求逆</b>
  求整数 `a_1, cdots, a_n` 模素数 `p` 的乘法逆元, 其中每个整数都与 `p` 互素.
  这是求逆元的一种离线算法. 记第 `i` 个<b>前缀积</b>为
  <span class="formula">
    `"pre"_i -= prod_(j=1)^i a_j (mod p)`.
  </span>
  注意若 `(a, p) = (b, p) = 1`, 则 `a^-1 b^-1 -= (a b)^-1 (mod p)`.
  于是前缀积的逆等于逆的前缀积, 即
  <span class="formula">
    `"pre"_i^-1 -= prod_(j=1)^i a_j^-1 (mod p)`.
  </span>
  利用这一性质设计算法如下:
  先使用 gcdExtended 计算所有数的乘积模 `p` 的逆元 `"pre"_n^-1`, 然后递推
  <span class="formula">
    `a_i^-1 -= "pre"_i^-1 * "pre"_(i-1) (mod p)`,<br>
    `"pre"_(i-1)^-1 -= "pre"_i^-1 * a_i (mod p)`.
  </span>
  时间复杂度为 `O(n + log p)`.
</p>
<div class="playground" value="{ a: [1, 1, 4, 5, 1, 4], p: 17 }">
<script type="text">
const { inv } = require('inv')
module.exports = function preProductInv (str) {
  const { a, p } = Playground.parse(str)
  // 初始化前缀积
  const pre = [1] // 0 个数的乘积是 1
  let last = 1
  a.forEach(x => {
    last *= x
    pre.push(last)
  })
  // 计算逆元
  let i = pre.length-1
  pre[i] = inv(pre[i], p)
  const ret = a.slice()
  // 注意数组 pre 的长度是 ret 的长度加 1
  for (i -= 1; i >= 0; --i) {
    ret[i] = pre[i+1] * pre[i] % p
    pre[i] = pre[i+1] * a[i] % p
  }
  return ret
}
</script>
</div>

<p class="algorithm">
  <b>试除法判断素数</b>
</p>

<div class="playground" value="{ n: 143 }">
<p>
143 = 11 * 13 不是素数
</p>
<script type="text" id="is-prime">
// 用试除法寻找 n 的最小素因子
// n = 1 时返回 1
function primeOf (n) {
  if (n % 2 === 0) return 2
  if (n % 3 === 0) return 3

  // 用 6n±1 试除
  const s = Math.sqrt(n)
  for (let p = 5; p <= s; p += 6) {
    if (n % p === 0) return p
    if (n % (p+2) === 0) return p+2
  }
  return n // n 是素数
}

// 判断 n 是否为素数
function isPrime (n) {
  return n > 1 && primeOf(n) === n
}

function demo (str) {
  const { n } = Playground.parse(str)
  return isPrime(n)
}

demo.primeOf = primeOf
demo.isPrime = isPrime
module.exports = demo
</script>
</div>

<p class="algorithm">
  <b>遍历所有因数</b>
</p>

<div class="playground" value="{ n: 30 }">
<p>30 的所有因数为 1, 2, 3, 5, 6, 10, 15, 30</p>
<script type="text">
module.exports = function factors (str) {
  let { n } = Playground.parse(str)
  const lo = []
  const hi = []
  const s = Math.sqrt(n)
  for (let d = 1; d <= s; ++d) {
    if (n % d === 0) {
      lo.push(d)
      hi.push(Math.floor(n/d))
    }
  }
  if (lo[lo.length-1] === hi[hi.length-1]) {
    lo.pop()
  }
  return lo.concat(hi.reverse())
}
</script>
</div>

<p class="algorithm">
  <b>因数分解 (暴力算法)</b>
  设 n 是正整数, 返回 n 的所有素因子及其次数.
</p>

<div class="playground" value="{ n: 252 }">
<p>252 = 2^2 * 3^2 * 7</p>
<script type="text">
const { isPrime } = require('is-prime')
module.exports = function factorNaive (str) {
  let { n } = Playground.parse(str)
  const d = {}
  for (let i = 2; i <= n; ++i) {
    while (n % i === 0) {
      d[i] = (d[i] || 0) + 1 // 空表项置为 1
      n = Math.floor(n / i)
    }
    // 可选, 起加速作用
    if (isPrime(n)) {
      d[n] = (d[n] || 0) + 1
      break
    }
  }
  return JSON.stringify(d)
}
</script>
</div>

<p class="remark">
  linux 系统可以用 <code>factor</code> 命令分解整数.
</p>

<p class="algorithm">
  <b>Pollard Rho因数分解, O(n^(1/4)), 1975</b>
  首先生成数列 `x_n`: 选取种子 `x_0`, `x_(n+1) = f(x_n) (mod n)`, 其中 `f(x) = x^2 + 1`. 设 `p | n`, 寻找整数 `x_i, x_j` 满足 `x_i = x_j (mod p)` 且 `x_i != x_j (mod n)`, 于是 `gcd(x_i - x_j, n)` 是 `n` 的非平凡因子.
</p>
<div class="playground" value="{ n: 65537 }">
<p>
65537 -&gt; prime<br>
9379 -&gt; 113<br>
32 -&gt; boom!
</p>
<script type="text">
const { gcd } = require('gcd')
// const { primeOf } = require('is-prime')
module.exports = function pollardRho (str) {
  let { n } = Playground.parse(str)
  const f = x => ((x * x % n) + 1) % n
  let [x, y] = [2, 1]
  while (x != y) {
    [x, y] = [f(x), f(y)]
    const g = gcd(n, Math.abs(x-y))
    if (1 < g && g < n) return g
    x = f(x)
  }
  return 'failed to factor, maybe prime?'
  // return primeOf(n) // fallback
}
</script>
</div>

<p class="algorithm">
  <b>Eratosthenes 筛法</b>
  求 n 以内的所有素数.
  注意到素数的倍数一定是合数,
  算法先将 bool 数组 <code>flag</code> 初始化为 false,
  一旦发现 j 是合数, 就令 <code>flag[j] = true</code>.
  时间复杂度为 `O(n log log n)`.
</p>

<div class="playground" value="{ n: 100 }">
<p>100 以内的素数为 2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97.</p>
<script type="text">
module.exports = function eratosthenesSieve (str) {
  let { n } = Playground.parse(str)
  const flag = []
  const primes = []
  for (let i = 2; i <= n; ++i) {
    if (flag[i]) continue
    primes.push(i)
    for (let j = i*i; j <= n; j += i) {
      flag[j] = true
    }
  }
  return primes
}
</script>
</div>

<table>
  <tr>
    <td>n</td>
    <td>`pi(10^n) = 10^n` 以内的素数个数</td>
  </tr>
  <tr>
    <td>7</td>
    <td>664,579</td>
  </tr>
  <tr>
    <td>8</td>
    <td>5,761,455</td>
  </tr>
  <tr>
    <td>9</td>
    <td>50,847,534</td>
  </tr>
</table>

<p class="algorithm">
  <b>Euler 线性筛</b>
  求 N 以内素数的一种线性时间的算法. 在 Eratosthenes 算法中,
  一个合数可能被反复判断: 如 12 被 2*6 排除, 也被 3*4 排除.
  现在规定每个合数都被自己的最小素因子排除, 从而 12 只能由 2 来排除,
  等等. 由于每个合数只被排除一次, 所以时间复杂度为 `O(N)`.
  <br/>
  算法使用的数组 <code>flag[N+1]</code> 同上.
  另有数组 <code>primes[M]</code> 用于保存已经得到的素数.
  可以参考上表确定合适的 M 值,
  如 <code>N = 1e9</code> 时可取 <code>M = 50847534</code>.
</p>

<div class="playground" value="{ n: 100 }">
<script type="text">
module.exports = function eulerSieve (str) {
  let { n } = Playground.parse(str)
  const flag = []
  const primes = []
  for (let i = 2; i <= n; ++i) {
    if (!flag[i]) primes.push(i)
    primes.some(p => {
      const k = i * p              // 一个合数
      if (k > n) return true       // k 太大了, 跳出循环
      // console.log(k, i, p)
      flag[k] = true               // p 是 k 的最小素因子, 将 k 标记为合数
      if (i % p === 0) return true // 详细的证明见下文
    })
  }
  return primes
}
</script>
</div>

<p class="proof">
  算法的关键一步:
  若素数 `p | i`, 则对任意素数 `q gt p`, `q` 都不可能是
  `k = i q` 的最小素因子, 所以算法内循环要在 <code>i % p == 0</code>
  时 break.<br/>
  现在证明这个算法确实使得每个合数都被自己的最小素因子排除.
  由于每个合数都可以写为自己的最小素因子和一个正整数之积,
  故只需证明对 `i = 2, 3, 4, cdots`, 任何形如 `k = i p`
  的合数都被 `p` 排除, 其中 `p` 是 `k` 的最小素因子.<br/>
  `i = 2` 时, 要使素数 `p` 是 `k = i p` 的最小素因子, 只有 `p = 2`.
  显然算法将 `4 = 2 xx 2` 排除, 因此 `i = 2` 时结论成立.<br/>
  `i gt 2` 时, 假设结论对所有小于 `i`, 大于等于 `2` 的正整数均成立,
  由于算法的外层和内层循环分别枚举整数 `i` 和素数 `p`, 所以合数
  `k = i p` 必然在某一时刻被排除.
  现在设 `p` 是 `k = i p` 的最小素因子, 但 `k = j q` 却被 `q` 所排除,
  其中 `q gt p` 是一素数. 设 `k_1 = j p`, 显然 `p` 也是 `k_1`
  的最小素因子.  注意到 `j lt i`, 由归纳假设, `k_1` 应该由 `p` 排除.
  然而由 `p | k = j q` 和 `p, q` 是不同素数知 `p | j`, 按算法,
  内循环应当执行 break 语句. 换言之, 不可能使内层循环继续迭代,
  从而排除合数 `k = j q`. 矛盾. 所以 `k` 必定由它的最小素因子排除.<br/>
  由数学归纳法, 结论对所有 `i = 2, 3, 4, cdots` 成立.
</p>

<p class="example">
  已知正整数 `n`, 求最大的正整数 `m | n`, 且 `m` 无平方因子,
  即不存在素数 `p` 使得 `p^2 | m`.
</p>

<p class="proof">
  为避免对 `n` 作因子分解, 可以求出 `2` 到 `sqrt n` 间所有素数的乘积
  `M`, 然后 `m = gcd(M, n)`.
  这里 `M` 是一个很大的数; 如果 `n` 的范围已知, 可以通过查表得到 `M`.
  如果 `M` 过大, 可以考虑每三个素数乘积求一次
  gcd, 再将结果相乘, 即
  <span class="formula">
    `m = gcd(2*3*5, n) * gcd(7*11*13, n) * cdots`
  </span>
</p>

<div class="p algorithm">
  <b>CRC 算法 (Cyclic Redundancy Check, 循环冗余校验)</b>
  是用于通信的校验算法.  该算法的输入是一个 `n+1` 位二进制串 `b` 和 一个
  `bbb F_2` 上的 `n` 次既约多项式 `f`.
  <br>
  所谓 `bbb F_2` 上的 `n` 次多项式是指, 系数为 0 或 1,
  系数上的加法和乘法是模 2 的 `n` 次多项式.
  `bbb F_2` 上所有次数不超过 `n` 的多项式共有 `2^n` 个, 它们构成一个域 `bbb
  F_(2, n)(x)`, 可以进行四则运算. 如乘法为
  <span class="formula">
    `(x^3 + x^2 + 1)(x^3 + x + 1)`
    `= x^6 + x^5 + x^4 + 3x^3 + x^2 + x + 1`
    `= x^6 + x^5 + x^4 + x^3 + x^2 + x + 1`.
  </span>
  以上过程简记为 1101 &times; 1011 = 111 1111. 又, `bbb F_(2, n)(x)`
  上的加减法恰好是异或运算:
  <span class="formula">
    1101 + 1011 = 1101 - 1011 = 0110.
  </span>
  所谓 `f` 是既约的, 是指 `f` 在 `bbb F_2` 上没有根, 这里就是说 0 和 1
  都不是 `f` 的根. 这容易办到, 只要保证 `f` 有奇数个系数为 1 的项,
  且常数项等于 1 即可.
  <br>
  当我们计算乘法所得的结果超过 `n` 次, 这个结果需要模 `f`, 使它落在 `bbb
  F_(2, n)(x)` 中. 比如, 当 `n = 4`, 取既约多项式 `f =
  x^4 + x + 1`, 则上面的结果化为
  <span class="formula">
    `x^6 + x^5 + x^4 + x^3 + x^2 + x + 1`
    `= (x^4+x+1)(x^2+x+1) + x^2 + x`
  </span>
  即, 111 1111 模 `f =` 1 0011 的最终结果是  `x^2 + x =` 110.
  这个求余数的过程很重要, 我们用竖式再算一遍 (注意减法就是异或运算,
  不需要借位):
<pre>
              111 --- 商 x^2 + x + 1 不重要, 可以不写
        ---------
1 0011 / 111 1111
         100 11
         -------
          11 001
          10 011
          -------
           1 0101
           1 0011
           ------
              110 --- 余数 x^2 + x
</pre>
  <br>
  现在可以介绍 CRC 的计算过程了:
  <span class="formula">
    `CRC(b, f) = x^n b(x) mod f(x) + x^n b(x)`.
  </span>
  将字节串 `b` 视作多项式, 首先将它乘以 `x^n`, 这相当于在串尾补 `n` 个 0.
  然后计算余数 `r(x) = x^n b(x) mod f(x)`, 这是一个次数不超过 `n` 的多项式.
  最后将 `x^n b(x)` 与 `r(x)` 相加, 即用 `r(x)` 替换串尾的 `n` 个 0 即可.
  <br>
  将数据 `b(x)` 和校验帧 `r(x)` 相连, 一起发送到接收方,
  接收方只需验证 `CRC(b, f) mod f(x) -= 0`, 就通过了 CRC 校验.
  若余数不为零, 则说明通信传输出现了问题.
</div>

<p class="remark">
  数字 `n` 称为 CRC 算法的位宽 (width),
  比如 CRC-8 的位宽等于 8, CRC-32 的位宽等于 32 等. 多项式 `f` 称为 CRC
  算法的生成多项式. 如 CRC-8 的生成多项式是 `x^8 + x^2 + x + 1`, CRC-32
  的生成多项式是
  <span class="formula">
    `x^32 + x^26 + x^23 + x^22 + x^16 + x^12 + x^11 + x^10 + x^8 + x^7 + x^5 + x^4 + x^2 + x + 1`.
  </span>
  `f` 用二进制表示时, 由于首项系数必为 1, 只记录 `n` 个低次项即可.
  如 CRC-8 和 CRC-32 的生成多项式分别记为 0x07 和 0x04C1 1DB7
</p>


下面这个计算 CRC-16 的例子来自 <a href="http://blog.csdn.net/liyuanbhu/article/details/7882789">CSDN</a>:
<pre>
#include &lt;stdio.h&gt;
#include &lt;stdint.h&gt;

const uint16_t POLY = 0x1021;

/* @param data 待校验的数据
 * @param len 数据包含的字节数
 * @param crc 初始的 CRC 余数
 */
uint16_t crc16(unsigned char *data, int len, uint16_t crc) {
  while (len--) {
    crc = crc ^ (*data++ &lt;&lt; 8); // fetch byte from memory, XOR into CRC top byte
    // loop for 8 bits
    for (int i = 0; i &lt; 8; i++) {
      if (crc &amp; 0x8000) // if b15 is set
        crc = (crc &lt;&lt; 1) ^ POLY;
      else // if b15 is off
        crc &lt;&lt;= 1;
    }
    crc &amp;= 0xffff; // ensure CRC remains 16-bit value
  }
  return crc;
}

int main() {
  uint16_t c1, c2;

  // 计算 data 的 CRC-16
  unsigned char data[] = {'1', '2', '3', '4', '5', '6', '7', '8', '9'};
  c1 = crc16(data, 9, 0xffff);
  printf("%04x\n", c1); // 29b1

  // 分两段计算 CRC-16
  c2 = crc16(data, 4, 0xffff);
  c2 = crc16(data + 4, 5, c2);
  printf("%04x\n", c2); // 29b1

  return 0;
}
</pre>

<p>
这里是 <a href="http://www.ip33.com/crc.html">CRC 在线计算</a>.
</p>

<ol class="algorithm" id="alg-cipolla">
  <b>Cipolla 算法: 求模 `p` 的平方根</b><br>
  [<a href="https://www.it610.com/article/1456933759239630848.htm">来自 olderciyuan</a>]
  设 `p` 为奇素数, `n` 为模 `p` 的二次剩余, 则 `n` 模 `p` 恰有两个平方根.
  Cipolla 算法能在 `O(log p)` 时间内求出这两个平方根.
  <li>取合适的 `a in bbb F_p`, 使得 `a^2 - n` 是二次非剩余.
    由于模 `p` 的二次非剩余恰有一半, 我们只要随机选取, 然后用 Euler 判别法,
    若 `(a^2-n)^((p-1)//2) = -1`, 即可判断它是二次非剩余.
  </li>
  <li>将 `a^2-n` 的一个平方根 `i = sqrt(a^2 - n)` 加到 `bbb F_p` 中,
    得到扩域 `bbb F_p(i)`. 这个域的特征为 `p`, 其中的元素形如 `A + Bi`,
    `A, B in bbb F_p`.
  </li>
  <li>可以证明 `(a+i)^(p+1) = n`, 因此 `n` 在扩域 `bbb F_p(i)`
    中有一对平方根 `+-(a+i)^((p+1)//2)`. 事实上, 这对平方根属于 `bbb F_p`,
    因此它们就是 `n` 模 `p` 的平方根.
  </li>
  <li>上述过程中应用快速幂算法, 总的时间复杂度为 `O(log p)`.
  </li>
</ol>

<ol class="proof">
  <li>先证 `i^p = -i`. 事实上
    <span class="formula">
      `i * i^(p-1)`
      `= i * (a^2-n)^((p-1)//2)`
      `= -i`.
    </span>
    最后一个等式成立是因为 `a^2-n` 为二次非剩余.
  </li>
  <li>Fermat 小定理: `AA a in bbb F_p`, `a^p = a`.</li>
  <li>由于 `bbb F_p(i)` 的特征为 `p`, 故 `AA a, b in bbb F_p(i)` 有
    `(a+b)^p = a^p + b^p`.
  </li>
  <li>现在可以算出:
    <span class="formula">
      `(a+i)^(p+1)`
      `= (a+i)(a+i)^p`
      `= (a+i)(a^p+i^p)`
      `= (a+i)(a-i)`
      `= a^2-i^2`
      `= n`.
    </span>
  </li>
  <li>下证 `n` 的平方根必定属于 `bbb F_p`. 设 `(A+B i)^2 = n`,
    即 `A^2 + B i^2 + 2A B i = n`. 比较两边虚部得 `A B = 0`, 因此 `A = 0`
    或 `B = 0`.<br>
    若 `A = 0`, 有 `B^2 i^2 = n`, 但 `B^2, n` 均为二次剩余, `i^2`
    为二次非剩余, 矛盾. 故 `B = 0`, 即 `A + B i in bbb F_p`.
  </li>
</ol>

<div class="playground" value="{ n: 8, p: 17 }">
<p>说明: `5^2 -= 8 (mod 17)`, 因此 8 模 17 的平方根为 `+-5`. 我们只显示其中一个根 (最小正根).</p>
<script type="text">
const { pow } = require('pow')

function legendre (n, p) {
  return pow(n, (p-1)/2, p) === 1 ? 1 : -1
}

// 在域 F_p(i) 上求 a = A+Bi 的 n 次方
function Field (p, ii) {
  this.p = p
  this.ii = ii
}

Field.prototype.mul = function (a, b) {
  return [
    (a[0] * b[0] + this.ii * a[1] * b[1]) % this.p,
    (a[1] * b[0] + b[1] * a[0]) % this.p
  ]
}

Field.prototype.pow = function (a, n) {
  if (n < 0) return NaN
  let ret = [1, 0]
  while (n) {
    if (n & 1) ret = this.mul(ret, a)
    a = this.mul(a, a)
    n >>= 1
  }
  return ret
}

// 求 n 模 p 的平方根
function cipolla (n, p) {
  n = (n % p + p) % p // 标准化
  if (n < 2) return n // sqrt(0) = 0, sqrt(1) = 1
  if (legendre(n, p) !== 1) return NaN

  // 寻找 a 使得 ii = a*a-n 为二次非剩余
  let a, ii, count = p
  do {
    --count
    a = 1 + Math.floor(Math.random() * (p-1)) // a = 1, 2, ..., p-1
    ii = (a*a - n + p) % p
    // console.log('a:', a, 'ii:', ii)
  } while (count && legendre(ii, p) === 1);
  if (!count) throw new Error('计算出错, 请检查输入')

  // 答案是 a+i 的 (p+1)/2 次方
  const res = new Field(p, ii).pow([a, 1], (p+1)/2)[0]
  return Math.min(res, p - res)
}

module.exports = function demo (str) {
  // 注: 若要列出模 p 的所有平方根, 不必使用 cipolla, 直接打表即可.
  const { n, p } = Playground.parse(str)
  return cipolla(n, p)
}
</script>
</div>

<script src="../../js/note.js?type=cs&loadmath=true"></script>
<script src="../../js/playground.js"></script>
</body>
</html>
